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I.  Introduction 

In  Chapter  2  a  state  space  formulation  of  the  Hoo  optimal  control  problem  is  given. 
Assuming  a  finite  interval  of  control,  the  problem  of  synthesizing  a  finite-interval  H^ 
controller  is  converted  into  an  optimization  problem  in  which  a  parameter  occurring  in  a 
boundary  value  problem  needs  to  be  maximized.  An  optimality  condition  for  the  maxi¬ 
mization  of  this  parameter  is  given.  The  proposed  method  makes  use  of  the  observer-based 
parametrization  of  all  stabilizing  controllers.  An  example  is  worked  out. 

An  important  problem  in  flight  control  and  flying  qualities  is  the  approximation  of  a 
complex  high  order  system  by  a  low  order  model.  In  Chapter  3,  for  a  given  reduced  order 
model,  we  define  the  correlation  measure  between  the  plant  and  the  model  outputs  to  be  the 
minimum  of  the  ratio  of  weighted  signal  energy  to  weighted  error  energy.  We  give  a  criterion 
for  the  evaluation  of  the  correlation  measure  in  terms  of  minimization  of  a  parameter 
occurring  in  a  two-point  boundary  value  problem.  Once  the  correlation  measure  for  a  given 
reduced  order  model  can  be  evaluated,  a  nonlinear  programming  algorithm  can  be  used  to 
select  a  model  which  maximizes  the  correlation  between  the  plant  and  model  outputs.  The 
correlation  index  used  can  be  regarded  as  an  extension  of  the  H oo  performance  criterion 
to  the  finite-interval  time-varying  case.  However,  the  usual  Hoo  problem  seeks  an  optimal 
controller,  whereas  our  problem  is  to  select  the  reduced  order  model  matrices  which  give 
the  best  correlation  index.  We  also  give  an  expression  for  the  variation  of  the  correlation 
owing  to  parameter  variations  and  pose  a  robust  model  reduction  problem.  The  utilization 
of  the  theory  is  demonstrated  by  means  of  some  examples.  In  particular,  a  problem  which 
involves  the  reduction  of  an  unstable  aircraft  model  with  structural  modes  is  worked  out. 
The  computations  for  Chapter  3  were  performed  by  Marc  Steinberg,  an  engineer  with  the 
Flight  Control  group. 
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II.  Synthesis  of  Finite-Interval  Hoo  Controllers 
by  State  Space  Methods 

l.  Introduction 

The  H0 o  optimal  control  theory  has  been  pioneered  by  Zames  [1]  and  important  con¬ 
tributions  have  been  made  by  FVancis  and  Doyle  [2,3].  Recent  work  [4]  indicates  that  the 
theory  has  important  applications  in  the  design  of  flight  control  systems. 

In  this  chapter  a  variant  of  the  problem  is  considered  in  terms  of  state  space 
formulation.  Optimization  routines  are  needed  for  the  synthesis  of  the  final  controller. 
The  formulation  is  based  on  considering  optimal  control  problems  with  finite  terminal 
time  in  which  the  cost  is  a  quotient  of  two  definite  integrals. 

Other  authors  have  considered  the  Hqo  problem  from  different  points  of  view.  In  [5] 
a  parametrization  of  all  stabilizing  controllers  that  achieve  a  specified  norm  bound 
is  given  in  a  specialized  case.  The  computation  of  the  controller  involves  the  solution  of 
two  Riccati  equations.  This  result  has  been  extended  to  the  general  case  in  [6].  In  [7]  the 
Hoo  problem  is  solved  by  introducing  a  generalized  algebraic  operation  called  conjugation. 
The  approach  again  yields  two  Riccati  equations  whose  solution  leads  to  the  synthesis  of 
a  controller.  In  [8]  a  certain  LQG  problem  with  a  side  constraint  on  the  #00 -norm  of 
the  closed  loop  transfer  function  is  solved.  In  this  approach  it  is  necessary  to  solve  three 
coupled  Riccati  equations.  In  special  cases  these  three  equations  can  be  reduced  to  two 
Riccati  equations. 

Our  approach  results  in  a  two-point  boundary  value  problem.  The  approach  has  the 
advantage  of  being  applicable  to  time- varying  systems  with  observer-based  controllers  and 
dynamic  controllers.  Ref.  9  contains  one  such  application  in  which  the  objective  is  to 
maximize  the  disturbance  rejection  capacity  of  a  time-varying  linear  system.  Also,  given 
a  controller  it  is  important  to  know  the  performance  measure  of  the  controller.  For  the 
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general  time-varying  system  with  a  given  controller,  the  parameter  A  of  Section  3  gives  a 
measure  of  the  performance  of  the  controller. 

Our  time-domain  approach  has  several  advantages  even  in  the  case  of  time-invariant 
systems.  First  of  all,  it  provides  an  alternate  new  approach  to  the  computation  of  finite- 
interval  Hoc  controllers.  The  ffoo  algorithms  usually  cannot  handle  time  domain  specifi¬ 
cations.  In  our  optimization  algorithm  it  is  possible  to  include  time  domain  constraints. 
Also  time  domain  approach  is  convenient  for  handling  parameter  uncertainties. 
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2.  State  Space  Formulation  of  the  Problem 
The  standard  i?oo  problem  can  be  stated  with  reference  to  Fig.  1  (p.  18).  In  Fig.  1 
w,u,z ,  and  y  denote  the  exogenous  input  (command  signals,  disturbances,  sensor  noises 
etc.),  the  control  input,  the  output  to  be  controlled,  and  the  measured  output,  respec¬ 
tively.  The  plant  G(s)  and  the  controller  K(s)  are  assumed  to  be  real-rational  and  proper. 
Partition  G  as 


The  equations  corresponding  to  Fig.  1  are 

z  =  G\\w  +  Gi2«,  V  =  G21W  4-  G22U,  u  =  Ky.  (2) 

The  standard  Hoo  problem  is  to  find  a  real-rational  proper  K  which  minimizes  the  H 
norm  of  the  transfer  matrix  from  to  to  2  under  the  constraint  that  K  stabilize  G. 

In  terms  of  state  space  equations  G(s)  is  written  as 

x  =  Ax  +  B\w  4-  B2U 

z  —  Ci  1  4-  DnW  ■+■  B12U  (3) 

y  ~  C2X  +  D21W  +  D22u- 

Doyle  [10]  showed  that  every  stabilization  procedure  can  be  realized  as  an  observer-based 
controller  by  adding  stable  dynamics  to  the  plant.  The  realization  of  the  observer- based 
controller  is  shown  in  Fig.  2  (p.  18)  where  the  stable  dynamics  added  is  represented  by 
Q(s),  with  Q(s )  proper  and  I  —  D22Q(oo )  invertible.  In  Fig.  2,  F  and  H  are  chosen  such 
that  A  +  B2F  and  A  +  HC2  are  stable.  Assume  that  Q(s)  is  described  by  the  minimal 
representation 

q  =  Aq  +  By,  u2  =  Cq  +  Dy.  (4) 


4 


NADC-90043-60 


Following  the  notation  of  [11],  define  the  following  quantities. 

/?i  =  -H  -  ( B2  +  HD22)(I  ~  DD22)~lD 
p2  =  B  +  BD22{I  -  DD22)~'b 

71  =  F  4-  (/  —  DD22)  b(C2  +  d22f ) 

72  =  -(I-DD22)~lC 

an  =  A  +  HC2  +  ( B2  4-  HD22)^\ 

(5) 

=  A  4-  B2F  —  Pi(C2  4-  D22F) 


012  =  (B2  +  H  1)22)72 


021  =  —&2(C2  4-  D22F ) 


022  —  A  —  BD22l2 

k  =  -{i-bD22)~lb. 

Then  the  closed  loop  system  is  given  by 

/  An  A12  A13  \  / x\  ( Bi  \ 

(  A21  A22  A23  1  I  ^  I  +  I  B2  I  w,  (6) 

\A3i  A32  A33  /  \q  J  \  ®3  J 


y  =  (I  —  D22k)  1[C2x  4-  D227ix  +  D22j2q  +  D2\w],  (7) 

z  =  C\X  4-  Dnw  +  £>12(71*  +  729  +  «y),  (8) 
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where 

Ajj  =  A  +  i?2«(J  —  -D 22k )  *^2 


A12  =  #27l  +  B2*(I  —  Di2*)  1-D227l 
A13  =  B272  +  B2*(I  ~  I?22«)  1B2272 
A21  =  /5l(J  —  £>22*)  ^2 
A22  =  <*n  4-  /?i(/  —  D22k)  1B>22l\ 

A-23  =  0112  +  fil(I  —  D22K)  1 D22I2 

(9) 

A31  =  02(1  ~  E>22k)  1C<2 

A32  =  <*21  +  02(1  ~  B22K)  1I?227l 

A33  =  <*22  +  ~  D22k)  1 D22I2 

Bi  =  Bi  -  B2D(I  -  D22K)  1D2 1 
B2  =  —(H  +  i?2-D).Z)2i 


B3  =  BD21. 

Consider  equations  (4)-(9).  Now  the  Hoc  control  problem  is  to  find  among  all  sets  of 
matrices  A,  B,  C,  and  D  which  give  a  stable  transfer  matrix  from  y  to  «2  (see  Fig.  2)  one 
for  which  the  Hoc-norm  of  the  transfer  matrix  from  w  to  z  is  minimized. 

The  above  problem  is  equivalent  to  the  following  problem.  Suppose  A  is  selected  to 
be  a  stable  matrix.  For  fixed  A ,  B,  C,  and  D,  let 

.  _  • ,  r 

where  the  superscript  *  denotes  matrix  or  vector  transpose.  Now  find  the  values  of  A,  B,  C, 
and  D  which  make  A  a  maximum.  The  initial  conditions  for  the  variables  x,x,  and  q  are 
of  course  zero. 

It  is  clear  that  the  Hoc-norm  of  the  transfer  function  from  w  to  z  is  l/y/X  and  the 
objective  is  to  minimize  the  Hoc -norm  by  choosing  a  controller. 
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The  input  w(t )  considered  in  the  above  problem  is  an  element  of  L 2(0,00).  However, 
in  many  physical  systems,  the  control  interval  is  finite.  For  example,  in  the  case  of  an 
advanced  fighter,  most  maneuvers  are  accomplished  in  the  course  of  a  few  seconds.  Thus, 
in  the  next  section  we  consider  an  approximate  Hoc  problem  in  the  sense  that  the  control 
interval  will  be  finite.  If  the  integration  limit  T  in,  say  equation  (13),  approaches  infinity, 
then  y/\  is  the  inverse  of  the  H oo-norm  of  the  transfer  matrix  from  w  to  z.  For  lack  of  a 
better  term,  we  call  this  a  finite-interval  Hoo  problem.  On  the  other  hand,  the  problem 
will  be  more  general  in  the  sense  that  time-varying  linear  systems  and  a  broader  class  of 
performance  indices  will  be  considered  in  Section  3. 

To  motivate  the  problem  considered  in  the  next  section,  let  x  =  (x*,x*,q*)* .  Equa¬ 
tions  (6)-(8)  are  written  as 

x  =  Ax  +  Bw,  x(0)  =  0,  w  =  u>,  (11) 

z  =  z  =  Cx  -f  Dw,  (12) 

where  the  matrices  A,  B,  C,  and  D  depend  on  A,  B,  C,  and  D.  Let  the  control  interval  be 
[0,  T].  For  fixed  A,  B,  C,  and  D  with  A  being  stable,  let 

Using  an  optimization  routine,  find  the  matrices  A,  B ,  C,  and  D  for  which  A  is  maximized. 

3.  Optimality  Conditions 

In  this  section  we  develop  conditions  for  determining  A  in  a  general  case  which  sub¬ 
sumes  the  problem  considered  at  the  end  of  Section  2.  These  conditions  will  be  developed 
for  time- varying  systems.  The  system  equations  are  given  by 

x  =  A(t)x  +  B(t)w,  x(t0)  =  0.  (14) 
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The  problem  on  hand  is  to  select  w  which  minimizes  the  performance  index  given  by 

T/_x  /^{Ix’RjX  +  x'R^  +  iw'Rjw}* 

«/( W )  —  mrr  (15) 

/to  U**WlX  +  x*W2w  +  ±w*W3w}  dt 

Note  that  the  performance  index  given  by  (13)  can  be  regarded  as  a  special  case  of  (15)  since 
z  =  Cx  +  Dw.  To  get  the  performance  index  of  (13),  set  Ri  =  R2  =  0,  Wj  =  C*C,  W2  = 
C*D,  and  W3  =  D*D  in  equation  (15).  In  (15)  we  assume  that  the  weighting  matrices 
Ri,R3,  Wx,  and  W3  are  symmetric  and  the  integrands  of  both  the  numerator  and  the 
denominator  are  nonnegative  for  each  w (<).  Further,  we  assume  that  there  is  some  w(t) 
for  which  the  denominator  is  positive.  Let  A  =  infw  J(w).  We  also  assume  that  R3  —  AW3 
is  nonsingular. 

Cost  functionals  of  the  form  of  (15)  are  the  subject  matter  of  this  report.  For  the  sake 
of  completeness,  we  derive  the  necessary  conditions  satisfied  by  an  optimal  w(<). 

Since  the  infimum  of  (15)  is  A,  we  have 
T 

f  {/x*R1x  +  x*R2w-|-  /w*R3w}dt 
Jt0  2  2 

T 

-a/  {/x’WjX+x’WjW+iw’Wawjdt  >0  (16) 

Jto  2  2 

for  all  (w,x)  which  satisfy  (14),  Thus,  if  w  minimizes  the  cost  functional  in  (15),  it  also 
minimizes  the  alternate  cost  functioned 
T 

Ji(w)=  /  {/x*(Rj -AW1)x  +  x*(R2-AW2)w  +  Iw*(R3_AW3)w}dt.  (17) 
Jt0  2  2 

The  necessary  conditions  for  optimal  w (<)  can  be  stated  as  follows. 

THEOREM  3.1.  Consider  the  system  given  by  (14)  with  the  performance  index  given  by 
(15).  Ifw(t)  minimizes  (15),  then  there  exists  an  adjoint  vector  if>(t)  such  that 

^  =  -A' V  +  (R,  -  AW,  )x  +  (R,  -  AW2)w,  t/>(T)  =  0,  (IS) 


8 


NADC-90043-60 


and 

w (t)  =  (R3  -  AW3)_1{B>  -  (R2  -  AW2)*x}.  (19) 

Proof.  To  give  a  short  proof,  consider  the  alternate  cost  functional  given  by  (17).  By 
the  maximal  principle  [12],  the  Hamiltonian  is  given  by 

w)  =  V»*(Ax  +  Bw)  -  {^x*(Ri  -  AWj)x 

+x*(R2  -  AW2)w  +  iw*(R3  -  AW3)w}.  (20) 


The  adjoint  vector  xp(t)  satisfies 


_  dH 
dt  dx 


(21) 


with  the  transversality  condition  t/i(T)  =  0.  Equation  (18)  is  obtained  from  (21).  Optimal 
w(t)  is  obtained  by  setting  dH /d w  =  0  and  is  given  by  (19).  □ 

Let  V,  =  R,  -  AW,  for  i  =  1, 2, 3.  We  have  a  two-point  boundary  value  problem  given 


by 


with 


xW  A-irv^v; 
rp)-  V  Vi  ~  V2Vr1V2 


BV^B* 

-A*  -  V2V^1B* 


x(<0)  =  0,  x!>{T)  =  0. 


(22) 


(23) 


We  now  show  that  the  minimum  value  of  (15)  is  the  least  positive  A  for  which  (22)-(23) 
has  a  solution  with  J^{Ax*WiX  +  x*W2w  +  |w*W3w}  dt  >  0. 

THEOREM  3.2.  Consider  the  boundary  value  problem  given  by  (22)  and  (23).  Let  A 
be  the  least  positive  value  for  which  the  boundary  value  problem  has  a  solution  with 
£{ix*WlX  +  x*W2w  +  jW*W3w}  dt  >  0,  where  w (t)  =  V3  ^B’V’  -  VJx}.  Then  A  is 
the  minimum  value  of  (15)  and  w  is  an  optimal  input. 

Proof.  From  Theorem  3.1  it  follows  that  if  w (t)  is  optimal,  then  the  boundary  value 
problem  (22)-(23)  is  satisfied  for  the  optimal  value  of  A.  Now  suppose  the  boundary  value 
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problem  is  satisfied  for  some  A  such  that  the  corresponding  solution  (x,  xp)  gives  the  de¬ 
nominator  of  (15)  a  positive  value  (with  w(<)  =  VJ'1{B*t/)  —  Vjx}).  We  show  that  the 
performance  index  corresponding  to  (x,  xp)  is  A. 

Let  (•,  •)  denote  the  standard  inner  product  in  a  real  Euclidean  space.  We  have 

(R3  -  AW3)w  =  BV  -  (R2  -  AW2)*x. 


Thus 


X*  T 

f  {(w,R3w)  -  A(w,W3w )}dt=  l  {(w,B»  -  (w,R£x)  +  A(w,W$x)}dt. 
Jto  Jto 

Since  Bw  =  x  —  Ax,  the  first  integral  on  the  right  side  of  (25)  is 
T  T 

j  (w,B*xp)dt=  [  {(x,V>)  -  (Ax,  xp))dt. 

Jto  Jto 

After  integrating  the  right  side  of  (26)  by  parts  and  utilizing  x(t0)  =  *P(T)  =  0, 

T  T 

[  (w,B *xp)dt  =  f  {(x,  Rix)  —  (x, R2 w)  +  A(x,  W3x)  +  A(x, W2w)}  dt. 
Jto  Jto 

Combining  equations  (25)  and  (27),  we  get 

f  {(x,  Rix)  +  2(x,  R2w)  +  (w,  R3w)}  dt  =  X  f  {(x,Wjx) 

Jto  Jto 

+2(x,  W2w)  +  (w,W3w)}  dt. 


(25) 


(26) 


(27) 


(28) 


Thus  A  is  the  cost  associated  with  (x,  xp).  Thus,  if  A  is  the  least  positive  value  for  which 
the  boundary  value  problem  (22)-(23)  has  a  solution  (x,  xp)  with  the  corresponding  denom¬ 
inator  of  (15)  being  positive,  then  x  must  be  an  optimal  trajectory.  □ 

If  the  system  and  weighting  matrices  axe  functions  of  a  finite  number  of  parameters, 
these  parameters  can  be  varied  to  maximize  A.  In  Section  2,  since  the  system  matrices  and 
the  weighting  matrices  depend  on  A,  B,C,  and  D,  an  optimization  routine  needs  to  be 
employed  with  respect  to  these  quantities  to  maximize  A. 


10 


NADC-90043-60 


4.  Optimality  Conditions  for  the  Maximization  of  A 
We  consider  again  the  time-invariant  H0 ©  problem.  In  this  section,  we  derive  a  con¬ 
dition  that  needs  to  be  satisfied  when  A  is  maximized.  For  this,  consider  equations  (14) 
and  (15).  Note  that  for  the  standard  problem  of  Section  2,  the  system  and  weighting 
matrices  depend  on  A,  B,  C,  and  D.  These  constitute  the  set  of  independent  variables. 
The  variations  in  the  system  and  weighting  matrices  can  be  explicitly  expressed  in  terms 
of  variations  in  A ,  B,  C ,  and  D.  However,  the  optimality  conditions  are  extremely  compli¬ 
cated  ;o  derive  in  such  a  case.  The  derivation  can  be  simplified  a  little  by  assuming  that 
£>22  =  0  (see  (3)).  However,  we  only  attempt  to  derive  the  basic  optimality  conditions 
here. 

Consider  equations  (22)  and  (23).  Let  A  =  A  —  =  BV^B*,  and  C  = 

Vi  —  "V^V^Vj.  Suppose  A,B,C ,  and  D  maximize  A.  Let  6A,6B,8C,  and  6D  denote 
elemental  perturbations  in  A,B,C,  and  D  respectively.  Also,  denote  the  corresponding 
perturbations  in  A,B,  C,x,rp,  and  A  by  6A,6B,6C,Xi,ipi,  and  fi  respectively.  Note  that  if 
A  is  a  maximum,  /x  =  0.  Thus,  we  have  the  following  set  of  equations. 


x  =  Ax  +  Bip, 

(29) 

ip  =  Cx-  A*ip, 

(30) 

eH. 

O 

II 

-e- 

v — ✓ 

II 

o 

(31) 

xi  =  Axj  +  Bip\  +  6Ax  +  SBip, 

(32) 

rpi  =  Cx  i  -  A*ipi  +SCx-  6A*ip, 

(33) 

Xi(<o)  =  tpi(T)  =  0. 

(34) 

F>om  (33),  we  have 

T  T 

f  X>1  dt=  /  {x*Cxi  -  x*A*ipi  +  x*SCx-  x*&4  ip}  dt.  (35) 

Jto  Jto 
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Also,  by  an  integration  by  parts 

rT  fT 

I  x*V>i  dt  =  -  /  {x*A*iJ>i+r{)*Bxf>i}dt. 

Jto  Jt0 

FVom  (35)  and  (36), 

fT  rT 

-  I  rp*Bxl>idt=  I  {x*Cxi  +  x*6Cx  -  x*6Atp)  dt. 

Jto  Jto 


Prom  equation  (30) 


Cx  =  t/>  +  A*V>- 


Note  that  C*  =  C.  Substituting  (38)  in  (37)  and  integrating  by  parts,  we  get 

fT  _  fT  ^  fT 

2  /  x*6ArJ>dt+  I  xfr'SBtfrdt  —  I  x*8C xdt  =  0. 

J  to  «/  to  •'to 


(36) 


(37) 


(38) 


(39) 


The  above  equation  needs  to  be  satisfied  for  all  elemental  perturbations  in  A,  B,  C, 
and  D. 


5.  A  Numerical  Example 

As  an  example  we  consider  the  tracking  problem  given  in  [2].  The  plant  is  given  by 


.  5  —  1 

P(s)  “  s(s-2)' 


(40) 


The  tracking  error  signal  is  r  —  v.  The  weighting  filter  W(s)  in  Fig.  3  (p.  18)  is  given  by 

3  +  1 


W(s)  = 


10s  + 1' 


(41) 


The  objective  in  [2]  was  to  choose  Ki(a)  and  Ki(s)  such  that  the  Hoc -norm  of  the  transfer 
function  from  w  to  v  is  minimized.  Our  objective  in  this  section  is  to  synthesize  u  using 
the  theory  of  this  chapter  such  that  the  minimum  of 


Jo10wa(0«ft 
fo°{(r-vf  +u2}dt 


(42) 
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is  maximized. 

Converting  the  plant  equations  to  state  space  form,  we  have 

X\  =  —  .lXj  +  w 


X2  =  « 


x3  =  2x3  +  u 


(43) 


r  =  .lit)  +  .09xi 


v  =  .5x2  +  .5x3. 


The  matrices  corresponding  to  equation  (3)  are  given  by 


The  matrices  F  and  H  are  chosen  such  that  A  +  B2F  and  A  +  HC2  are  stable.  The 


choice  is  the  same  as  that  in  [2]  and  is  given  by 

(0  0\ 

F  =  (  0  .5  -4.5),  H=  0  1  . 

\0  -9 ) 

Assume  that  Q(s)  is  described  by  the  three  dimensional  system 

q  =  Aq  +  By, 


u2  =  Cq  +  Dy. 


(44) 


Let  x  =  ( xi  X2  X3  )*  .  Then  the  state  equations  for  the  finite-interval  Hoo  problem 
become 


x  =  ( A  —  B2DC2)x  +  B3{F  +  DC2)x  —  S2C9  +  (^i  —  B2DD2i)w,  (45) 

i  =  {A  +  hc2  +  b2f  +  B2DC2)x  -  ( HC2  +  B2DC2)x 


—B2Cq  —  (HD2i  +  B2DD2\)w, 


q  =  Aq  4-  BC2{x  -  x)  +  BD21W, 
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with  the  initial  conditions  being  zero.  The  performance  index  is 

_ fo10^dt _ 

-10  f  (.lu>  +  .09xi  -  .5x2  -  -5x3 )2 

1  +[Fx  -  (Cq  +  DC2x  -  DC2x  +  DDiiw)}'2 


jcft 


Assuming  values  for  A,  B,  C,  D,  we  can  find  A  using  the  theory  given  in  Section  3. 

Let  $  =  [  f*11  ^12  )  be  the  transition  matrix  corresponding  to  (22).  Satisfaction  of  (23) 

\$21  $22  J 

gives  rise  to  the  condition  that  det($22(10))  =  0.  Thus  A  is  found  by  making  use  of  a  sign 
change  of  det($22(10))  over  a  range  of  values  of  A.  In  our  numerical  experiments,  much  of 
the  computer  execution  time  was  consumed  by  the  calculation  of  A  for  a  given  controller. 
Efforts  are  under  way  to  make  the  computation  of  A  more  efficient. 

The  transition  matrix  $(10)  was  found  in  this  case  using  the  following  formula  [13]. 
Let  h  =  10/28.  Represent  the  system  matrix  in  (22)  by  M.  Then 

2s 

=  +  +  +  .  (49) 


Using  the  above  procedure,  we  can  iterate  on  A,  B,  C,  and  D  to  maximize  A.  Note  that 

28 

once  $(h)  is  calculated,  only  eight  repeated  squarings  are  needed  to  evaluate  {$(/i)}  . 

Initially  the  following  values  were  assumed  for  the  control  matrices: 

.  /“2  0  0  \  /  1  -i\ 

A  =  [  0  -2  0  ,B=  -1  1  ,C  =  (1  1  1 ) ,  £>  =  ( 1  1). 

\  0  0  -2  7  \  1  -1/ 


Using  the  Rosenbrock  hill  climbing  algorithm  [14],  the  elements  of  the  matrices  were  varied 
to  maximize  A.  The  algorithm  usually  leads  to  only  local  maxima.  Note  that  Q(s)  is  stable 
if  and  only  if  A  is  stable.  This  was  not  introduced  as  a  constraint  in  the  optimization 
algorithm  since  the  unconstrained  run  yielded  a  stable  A.  The  Fortran  program  was  run 
on  a  Zenith  Z-248  personal  computer  in  double  precision  using  the  Microsoft  Optimizing 
Compiler  Version  4.01.  A  local  maximum  of  A  =  14.8  was  obtained  for  the  following  values 
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of  A,B,C ,  and  D: 


/  —2.04 

.318 

.023  \ 

(  .945 

9.973  \ 

A  =  -.026 

-1.632 

-.028  , 

B  =  [  -1.046 

33.028 

V  -.052 

.358 

-2.054  J 

\  .946 

-1.056  ) 

C  =  ( .944  1.48  1.018),  X>  =  (.986  41.92). 


After  several  runs  with  various  initial  values  for  A,  B,  C,  and  D,  the  value  of  Amax  =  14.8 
could  not  be  bettered. 


The  two  compoents  of  Q(s)  are  given  by 


QM 

Qi(*) 


.986(s  +  3.26)(s  +  2.03)(s  +  .75) 
(s  +  1.68)(s2  +  4.05s  +  4.1)  ’ 

41.92(s  +  2.04)(s2  +  5.04s  +  6.58) 
(s  +  1.68)(s2  +  4.05s  +  4.1) 


(50) 


It  was  reported  in  [2]  that  Q2(5)  is  unconstrained  and  may  be  taken  as  zero.  To 
simulate  this  condition,  we  set  the  second  columns  of  the  optimal  B  and  D  equal  to  zero. 
The  first  positive  value  of  A  for  which  det($22(10))  changed  sign  in  this  case  was  still 
observed  to  be  14.8. 

A  few  comments  on  txie  numerical  method  axe  in  order.  Since  the  computation  of  A 
consumes  most  of  the  execution  time,  further  research  needs  to  be  done  to  find  an  alternate 
method  to  evaluate  A  more  accurately  and  efficiently.  Also,  the  value  of  A  is  evaluated 
in  the  above  example  by  starting  with  an  initial  value  and  incrementing  it  in  steps  of  0.2 
until  a  change  in  the  sign  of  the  determinant  is  observed.  Thus  the  exact  value  of  A  differs 
from  the  computed  value  by  at  most  0.2.  This  sort  of  inaccurate  evaluation  of  A  may 
prematurely  terminate  the  optimization  routine  which  seeks  to  maximize  A. 


6.  Conclusions 

A  design  methodology  for  the  synthesis  of  finite-interval  Hoo  controllers  is  presented 
using  state-space  methods.  Using  observer-based  controller  parametrization,  an  optimiza¬ 
tion  problem  is  formulated.  A  measure  of  performance  for  a  given  controller  is  defined 
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in  terms  of  the  least  value  of  a  parameter  occurring  in  &  two-point  boundary  value  prob¬ 
lem.  Optimality  conditions  for  finding  the  measure  of  performance  for  a  given  controller 
are  given.  The  optimization  problem  seeks  to  maximize  the  measure  of  performance.  An 
example  is  given. 

Note:  This  chapter  is  based  on  a  paper  which  will  appear  in  the  AIAA  Journal  of 
Guidance,  Control,  and  Dynamics  under  the  same  title. 
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Fig.  1  The  Standard  Block  Diagram 


Fig.  2  The  Observer-based  Controller  Parametrization 


Fig.  3  Block  Diagram  for  the  Tracking  Problem 
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III.  Model  Reduction  with  a  Finite-Interval 
Hoc  Criterion 

1.  Introduction 

Model  reduction  is  an  important  problem  in  the  case  of  airplanes  with  significant 
aeroservoelastic  dynamics.  The  original  model  in  such  cases  is  of  high  order  and  thus,  the 
resulting  controller  will  have  a  complex  structure,  especially  if  it  uses  full  state  feedback. 
Also  for  highly  augmented  aircraft  with  flight  and  propulsive  controls,  it  is  useful  to  develop 
low  order  models  to  analyze  flying  qualities. 

If  the  aim  is  to  design  a  low  order  controller  for  a  high  order  plant,  there  are  at  least 
three  broad  approaches  to  achieve  this.  A  general  account  of  these  three  approaches  is 
given  in  [1],  The  so  called  direct  design  methods  assume  a  stabilizing  controller  of  fixed 
degree  and  seek  to  find  the  controller  that  maximizes  a  quadratic  performance  index  (  see 

12,3)). 

Another  approach  is  to  get  a  high  order  controller  by  some  design  technique,  such  as 
LQG  or  tfoo,  and  then  to  approximate  the  high  order  controller  by  a  low  order  one  which 
possesses  certain  desirable  properties.  This  approach  is  the  subject  matter  of  [1]  and  the 
pertaining  literature  is  referenced  in  that  paper. 

The  third  approach  is  to  approximate  the  high  order  plant  by  a  low  order  one.  Then 
a  low  order  controller  is  designed  and  used  to  control  the  original  plant.  In  this  chapter 
we  concentrate  on  this  approach  and  consider  the  problem  of  approximating  the  original 
plant  by  a  low  order  model  in  an  optimal  sense.  This  problem  has  been  treated  recently 
by  several  researchers  under  a  variety  of  approximation  criteria  and  we  refer  the  reader 
to  [4]  for  the  relevant  references.  Although  no  computational  results  are  given,  [4]  gives  a 
sufficient  condition  which  characterizes  reduced  order  models  satisfying  an  optimized  L2 
bound  as  well  as  a  prespecified  H bound.  The  reduced  order  model  is  expressed  in  terms 
of  solutions  of  four  coupled  algebraic  Riccati  equations. 
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We  now  state  the  main  problem.  For  the  sake  of  generality,  we  pose  it  for  time- varying 
systems.  Let  the  plant  be  described  by 

xp  =  Ap(t)xp  +  Bp(i)u,  xp(t  o)  =  0,  (1) 

yp  =  Cp(t)xp  +  Dp(t)u ,  (2) 

where  xp(t),u(t),  and  yp(t)  denote  the  plant  state  vector,  the  control  vector,  and  the  plant 
output  vector  respectively.  Let  the  reduced  order  model  which  approximates  the  plant  be 
chosen  to  be 

im  =  Am(t)xm  +  Bm(t)u,  Xm(to)  =  0,  (3) 

Vm  =  cm{t)xm  +  Dm(t)u,  (4) 

where  xm(t)  and  ym(t)  denote  respectively  the  state  vector  and  the  output  vector  of  the 
reduced  order  model. 

For  given  Am(t),Bm(t),Cm(t),  and  Dm(t ),  let  u  be  chosen  such  that  the  correlation 
index  given  by 

T  (  0  J 

Ito  2 (yp  -  ymTQ(t){yp  -  ym)dt 

is  minimized.  The  superscript  *  denotes  matrix  or  vector  transpose.  Let  this  minimum 
value  be  denoted  by  A.  Thus  u  represents  the  worst  input  and  A  gives  a  measure  of  the 
worst-case  correlation  between  the  plant  output  and  the  model  output.  The  problem  is  to 
choose  Am(t),Bm(t),  Cm(<),  and  Dm(t)  such  that  A  is  maximized. 

Since  (5)  represents  the  ratio  of  weighted  signal  energy  to  weighted  error  energy,  the 
above  problem  may  be  regarded  as  a  modified  problem  except  for  a  few  differences. 
We  consider  time-varying  systems  and  in  our  case  the  interval  of  control  is  finite.  There 
are  extensions  of  the  results  to  the  finite-interval  time- varying  case  [4].  However,  our 
approach  is  different  and  is  based  on  considering  the  inherent  two-point  boundary  value 
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problem.  Also,  the  general  aim  of  Hoo  problems  is  the  design  of  an  optimal  controller, 
whereas  in  this  chapter  we  are  interested  in  the  selection  of  model  matrices.  It  is  necessary 
in  our  case  to  use  nonlinear  programming  algorithms  in  order  to  select  the  model  matrices 
which  maximize  A.  In  [5-7],  we  derived  some  results  which  aid  in  the  selection  of  a  controller 
which  maximizes  the  worst-case  performance.  The  results  of  [5]  are  presented  in  Chapter 
2  and  these  will  be  utilized  in  Section  2  of  this  chapter. 

In  the  case  of  time-invariant  systems,  a  nonlinear  programming  algorithm  can  be 
used  to  find  at  least  a  local  maximum  of  A.  For  the  time-varying  case,  the  matrices 
Cm(t),  and  Dm(t)  need  to  be  expressed  in  terms  of  basis  functions  and  a 
nonlinear  programming  algorithm  needs  to  be  used  to  maximize  A  with  respect  to  the 
coefficients  of  the  basis  functions. 

We  do  not  require  the  plant  and  the  model  to  be  open  loop  stable.  This  is  significant 
since  many  of  the  modern  aircraft  have  open  loop  unstable  poles.  We  show  in  Section  4  by 
means  of  examples  that  the  method  is  indeed  applicable  to  such  cases.  There  is  yet  another 
advantage  of  our  method.  One  of  the  criticisms  in  the  approach  of  getting  a  low  order 
model  from  a  high  order  plant  is  that  the  satisfactory  approximation  of  the  plant  requires 
some  knowledge  in  advance  of  the  controller  [1].  Since  we  maximize  the  correlation  between 
the  plant  and  model  outputs  for  the  worst  possible  input,  the  correlation  in  the  case  of 
any  other  controller  is  bound  to  be  better.  Thus,  our  method  furnishes  a  satisfactory 
approximation  without  requiring  an  a  priori  knowledge  of  the  controller. 

We  now  give  a  summary  of  the  results  of  the  chapter.  In  Section  2,  conditions  that 
characterize  the  worst  input  are  derived  for  a  given  model.  A  two-point  boundary  value 
problem  needs  to  be  solved  for  the  least  positive  A  to  obtain  the  worst-case  correlation 
between  the  outputs  of  the  plant  and  the  model.  A  nonlinear  programming  algorithm  can 
then  be  used  to  find  the  model  matrices  which  maximize  A. 

The  stability  and  control  derivatives  of  aircraft  are  subject  to  variations  and  it  is  also 
not  possible  to  determine  these  exactly  from  wind  tunnel  data.  There  is  already  some 
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interest  in  robust  model  reduction  techniques  [8].  In  Section  3  we  formulate  a  robust 
model  reduction  problem  and  derive  an  expression  for  the  variation  of  correlation  between 
the  plant  and  model  outputs  as  a  functional  of  the  variations  in  system  parameters.  This 
value  gives  an  idea  of  the  robustness  of  the  approximate  model  and  can  aid  in  the  choice 
of  a  reduced  order  model  with  a  specified  level  of  robustness. 

In  Section  4  some  examples  are  worked  out  and  details  about  the  computational 
algorithm  utilized  are  given.  Correlation  between  the  plant  and  the  model  is  shown  via 
time  and  frequency  response  plots.  In  order  to  keep  the  examples  as  simple  as  possible, 
we  do  not  consider  the  robust  model  reduction  problem  in  the  case  of  these  examples. 

Finally,  certain  conclusions  are  given  in  Section  5. 

2.  Computation  of  A  for  a  Given  Reduced  Order  Model 

Assume  that  the  matrices  Am(t),  Bm(t),  Cm(t),  and  Dm{t)  are  given.  In  this  section, 
we  characterize  A  as  the  minimum  positive  value  for  which  a  certain  two-point  bound¬ 
ary  value  problem  has  a  nontrivial  solution.  Also,  we  derive  a  computationally  useful 
characterization  of  A. 

Letting 

x*  =  (x'P  O* > 
y  =  yP-  ymj 

^>=(o  l)- 

*«-(&)• 

c(t)  =  (  Cp  -Cm), 

and 

D(t)  =  Dp  —  Dm, 


(6) 

(7) 

(8) 

(9) 

(10) 

(11) 
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we  can  write  (l)-(4)  as 

x  =  A(t)x  +  B(t)u ,  i(<o)  =  0, 

y  =  C(t)x  +  D(t)u. 

The  correlation  index  given  by  (5)  can  be  put  in  the  form 

“  f^{jX’W,x  +  x’W2u  +  \u"W,u}dt' 

The  problem  is  to  characterize  u(t)  that  minimizes  (14). 

Let  A  =  infu  J(u).  We  assume  that  for  all  t,  R(t)  —  A W3(t)  is  invertible,  and 


R(t)  >  0, 


fW1(t)  W2(t) 
\W2-(t)  W3(t) 


^  >0. 


Equations  (15)  and  (16)  guarantee  that  the  numerator  and  denominator  of  (14)  are  non¬ 
negative  for  any  u.  The  necessary  conditions  that  characterize  the  worst  input  can  be 
stated  as  follows. 

THEOREM  2.1.  Consider  the  system  given  by  (12)-(14).  If  u(t )  minimizes  the  correlation 
index  given  by  (14),  then  there  exists  an  adjoint  vector  tft(t),  not  identically  zero,  such  that 


=  -A* ip  -  AWix  -  AW2u,  xJ>(T)  =  0, 
dt 


u{t)  =  (R-  A +  AW2*x). 
Proof.  For  a  proof,  see  Theorem  3.1  of  Chapter  2. 


A-A  +  A B(R  -  AW3)_1  W2*, 


B  =  B{R-XW3)-lB\ 
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c  =  -\Wi  -  a 2w2(r  -  a w3)~lw;. 


(21) 


Thus,  we  have  a  two- point  boundary  value  problem  given  by 


(22) 


with 


x(t0)  =  0,  rP(T)  =  0.  (23) 

The  following  theorem  follows  from  Theorem  3.2  of  Chapter  2. 

THEOREM  2.2.  Let  (x,V>)  satisfy  the  boundary  value  problem  given  by  (22)  and  (23) 
for  the  least  positive  A  such  that  fto{%x*W\X  +  x’W^u  +  iu’Wau}  dt  >  0,  where  u  = 
(R  —  A W 3 )—1  {B*xp  +  \W2x}.  Then  A  is  the  minimum  value  of  the  index  given  by  (14)  and 
u  is  the  worst  input. 


In  [5]  and  [6],  a  computational  technique  which  utilizes  the  transition  matrix  associated 
with  (22)  is  given.  This  technique  is  also  presented  in  Chapters  2  of  this  report.  In  this 
chapter  we  use  an  alternate  technique.  This  technique  is  more  stable  numerically.  The 
theory  behind  the  technique  is  given  below. 

Let  $(t,r)  be  the  transition  matrix  associated  with  (22).  Then  we  have 


(24) 


(25) 


(26) 


(27) 
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Since  x(to)  =  ip(T)  =  0, 

Cn  x(T)  =  ui2tl>(to)y 
C21  x(T)  =  V22lJ>(to)- 

Since  the  equations  in  (28)  are  linearly  dependent,  (28) 
V>(<o)  and  x(T)  if  and  only  if 


(28) 

has  a  nontrivial  solution  for 


det 


(<ll  •'12)=0. 

\  C21  Vw  ) 


(29) 


Thus,  we  can  characterize  A  as  the  least  positive  value  for  which  (29)  holds. 

We  can  determine  A  by  doing  a  search  over  a  range  of  positive  values  and  picking  the 
first  value  at  which  the  determinant  in  (29)  changes  sign.  We  give  more  details  on  this  in 
Section  4. 


3.  Robust  Model  Reduction 

In  this  section  we  formulate  a  robust  model  reduction  problem.  The  aim  is  to  choose 
the  best  reduced  order  model  under  parameter  variations.  We  derive  an  expression  for 
the  variation  in  the  correlation  measure  A  in  terms  of  variations  in  the  system  matrices. 
For  simplicity  of  analysis,  we  assume  that  Dp(t)  and  Dm(t)  in  (11)  are  zero,  which  makes 
D(t)  =  0. 

Consider  (1)-(10).  The  system  equations  are  given  by 

x  =  A(t)x  +  B(t)u ,  x(<o)  =  0,  (30) 

y  =  C(t)x.  (31) 

We  can  write  (5)  as 

/I  hu'(t)R(t)u(t)dt 

Si  \z‘it)C'(tW)C(t)x{t)it. 

For  given  Am(t),Bm(t),  and  Cm(t),  let  A  be  the  minimum  of  the  correlation  index 
in  (32)  over  u(t).  Let  the  elemental  variations  in  Ap(t),  Bp(t),  and  Cp(t)  be  denoted  by 
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6Ap(t),6Bp(t ),  and  SCp(t)  respectively.  Let  &4(t),ffi(t),  and  SC(t )  be  the  variations  in  the 
matrices  A(t),B(t),  and  C(t)  corresponding  to  the  elemental  variations  6Ap(t),dBp(t),  and 
6Cp(t).  Notice  that 


m  =  (SAf)  ®). 

®w=(T)’ 


6C(t)  =(K7,(t)  0). 


(35) 


Let  n  denote  the  variation  in  A  caused  by  &4,  SB,  and  SC.  Now  the  robust  model 
reduction  problem  can  be  stated  as  follows. 

Robust  model  reduction  problem.  Find  Am(t),  Bm(i),  and  Cm(t)  such  that 

JtT  \u*Rudt 

inf  - 

“  Si  \x*C*QCxdt 
is  maximized  with  the  side  constraint 


l/*A|  <  Ho  for  all  ||&4(t)||  <  a(t),  ||6B(t)||  <  b(t),  and  ||«7(<)||  <  c(t),  (37) 

where  a(t),  b(t),  and  c(t)  are  suitably  chosen. 

We  now  derive  an  expression  for  /i  in  terms  of  SA(t),£B(t),  and  SC(t).  For  given 
Am(t),  Bm(t),  and  Cm(t),  let  u  minimize  the  index  in  (32).  In  the  following,  we  suppress 
the  dependence  of  the  matrices  on  t  for  simplicity  of  notation.  From  (19)-(23),  we  get  the 
following  boundary  value  problem  which  needs  to  be  satisfied  by  the  corresponding  pair 
(x,rp). 

i=Ax  + 

=  -A C'QCx  - 


26 


(38) 

(39) 


NADC-90043-60 


x(t0)  =  xl>(T)  =  0.  (40) 

Let  xi  and  xpi  represent  the  variations  in  x  and  tp  due  to  SA,SB ,  and  SC.  From  (38)- 
(40),  we  have  the  following  equations  satisfied  by  xx  and  tp\. 

ii  =  Axx  +  &A  x  +  BR~lB'r)> j  +  ( BR-'SB *  +  SB  R~lB*)tP,  (41) 

tpi  =  -nC'QCx  -  A (6CmQC  +  C*QSC)x  -  \C*QCxx  -  A'fa  -  &4>,  (42) 


xi(to)  =  *Pi(T)  =  0. 


Theorem  3.1.  Consider  (38)-(43).  Then  the  variation  in  A  is  given  by 

-  /tT  xp'SA  x  dt  —  /tT  tp*B*R~lSB  ip  dt  —  A  JT  x*C*Q«7  x  d< 


fi  =  — 


Sl\x*C*QCxdt 


Proof.  From  (42),  we  get 


[T  x*xj>1dt  = -fi  f x*C*QCxdt-X  f x*(SC*QC  +  C*QSC)xdt 

Jto  j  to  Jto 

-A  f  x*C*QCxi  dt  -  f  x*A*xpi  dt  —  f  x'SA'xpdt. 

J  tQ  "to  "to 


Also,  by  an  integration  by  parts  and  by  (38),  (40),  and  (43), 


Since 


tT  fT  tT 

/  xipi  dt  —  —  /  xmAxPxdt-  /  rp*BR-1  B*xpi  dt. 
Jto  Jto  Jto 

[  x'(SCmQC  +  C*QSC)xdt  =  2  f x'C'QSC  xdt, 
Jto  Jto 


from  (45)  and  (46),  we  get 


i  /  x*C*QCxdt  +  2A  [T x'C'QSC  xdt 
Jto  Jto 

+  A  /  x*C*QCxi  dt  +  [T x*6A*iPdt=  C tP*BR-lB*rpxdt. 

J  to  "t$  "to 
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FVom  (39), 

tT  tT  •  * 

A  /  x'C'QCx i  dt  =  -  (4>  +  A*xl>)  Xl  dt.  (49) 

Jto  Jto 

Integrating  the  first  term  of  the  integrand  from  the  right  side  of  (49)  by  parts,  and  using 
(40)  and  (43),  we  get 

A  [T  x*C*QCx1dt=  f  xJ>*6Axdt+  f  ^BIC'B'xhdt 
Jto  Jto  Jto 

rT 

+  /  rl>*{BR~l8B' +  6B  R~'B*)xl>dt.  (50) 

Jto 

Incorporating  (50)  in  (48)  and  using  the  fact  that 

rjy 

f  \})*(BR~X5B*  +  8B  R~lB*)tl)dt  =  2  f  tp*BR~x6B  dt,  (51) 

Jto  Jto 

we  get  (44).  □ 

Using  (44),  the  variation  in  the  correlation  measure  A  owing  to  parameter  variations 
can  be  computed  for  any  given  Am(t),  Bm(t),  and  Cm(t). 


4.  Numerical  Examples 

In  this  section  we  consider  only  time-invariant  examples.  The  systems  in  the  examples 
can  be  put  in  the  form 


x  =  Ax  +  Bu ,  x(0)  =  0,  x*  —  ( z!  x^  )*  , 


(52) 


with  the  correlation  index 


JM  =  l  iu'Rudt 

JW  tT  IrT't 


(53) 

So1  \x*C*QCx  dt 

where  A,B,  and  C  are  given  by  (8)-(10).  For  given  Am,Bm ,  and  Cm,  let  A  =  inftt  J(u). 
To  recap  the  procedure  for  finding  A,  let 


/  A  BR~lB*  \ 
_^-AC*QC  -A*  )' 


(54) 
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and 

“»<?)' -ft  <55> 

<;:)•  (56> 

From  (29),  A  is  given  by  the  first  positive  value  for  which 


Now  we  iterate  on  the  matrices  Am,Bm,  and  Cm  using  a  nonlinear  programming  algorithm 
to  maximize  A. 

There  are  two  primary  computational  algorithms  that  are  needed  to  use  this  method 
of  model  reduction.  They  are  a  nonlinear  global  optimization  routine  to  maximize  A  and 
a  relatively  fast  routine  to  compute  A.  Currently  A  is  determined  by  finding  the  smallest 
positive  value  for  which  (57)  holds.  Due  to  the  oscillatory  nature  of  the  value  of  the 
determinant  as  a  function  of  A,  for  suitable  weighting  matrices  in  (53),  A  was  originally 
calculated  starting  with  an  initial  value  of  A  =  0.1  and  incrementing  by  0.2  until  the 
determinant  in  (57)  changed  sign.  While  this  method  yielded  accurate  results,  it  also  used 
excessive  amounts  of  computational  time. 

To  speed  this  process  up,  two  modifications  were  made.  First,  A  was  incremented  by 
large  steps  until  the  value  of  the  determinant  was  less  than  a  percentage  of  its  initial  value. 
At  this  point,  small  increments  in  A  were  applied  until  the  determinant  changed  sign.  This 
technique  was  successful  in  this  case  because  the  absolute  value  of  the  determinant  never 
increased  beyond  a  very  small  fraction  of  its  initial  value  after  the  first  zero  crossing.  The 
second  modification  was  adaptively  scaling  down  the  input  weighting  matrix  R  so  that  the 
values  of  A  were  consistently  in  the  range  of  10  to  20.  This  modification  gives  large  enough 
values  for  accuracy  and  small  enough  values  to  decrease  computational  time. 

The  second  necessary  algorithm  is  a  nonlinear  global  optimization  routine.  We  are 
still  in  the  process  of  developing  an  algorithm  which  satisfies  both  speed  and  accuracy 
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requirements.  For  now,  however,  two  methods  were  used  to  test  our  theory.  For  Example 
1,  we  used  a  deterministic  tunneling  technique  [9].  This  technique  in  our  case  starts  with  a 
modified  version  of  the  Rosenbrock  constrained  hill  climbing  algorithm  [10],  searches  for  a 
better  point  than  the  current  local  maximum,  and  then  restarts  the  hill  climbing  algorithm 
from  there.  While  this  method  did  converge  to  the  global  maximum,  it  also  required  an 
excessive  number  of  iterations.  The  method  used  for  Example  2  was  a  multi-start  hill 
climbing  algorithm  with  the  starting  point  chosen  by  truncating  the  original  system  as 
well  as  by  other  model  reduction  techniques.  This  method  will  in  general  not  converge 
to  the  global  maximum  without  an  excellent  starting  point,  but  it  will  find  a  good  local 
maximum  for  a  decent  starting  point. 

The  Rosenbrock  hill  climbing  algorithm  and  the  algorithm  to  compute  X  were  written 
in  PC-MATLAB  and  run  on  a  Zenith  Z-248  personal  computer. 

Example  1.  Simple  illustrative  examples 

We  will  first  show  the  reduction  of  a  stable  second  order  system  and  an  unstable 
second  order  system.  The  system  is  of  the  form 


xp  =  ApXp  +  Bpu,  xp(  0)  =  0, 
yp  —  CpXp , 


where  At 


=  f°r  4 he  stable  case,  and  Ap  = 

case.The  other  matrices  are  Bp  =  ^  ^  and  Cp  =  ( 1  0  ) . 
Our  first  order  model  equations  are 


(i°o  -9) 


(58) 

(59) 

for  the  unstable 


%m  —  amxm  "l"  ^tn(O)  —  0,  J/m  —  CmXm .  (60) 

The  weighting  matrices  in  (53)  are  R  =  .01,  Q  =  1,  and  the  final  time  T  is  taken  to  be  2 
seconds.  Our  optimization  technique  is  the  aforementioned  tunneling  algorithm. 
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In  the  stable  case,  the  initial  values  were  am  =  —1,  bm  =  1,  and  cm  =  1.  The  value  of  A 
increased  from  an  initial  value  of  3.4  to  a  maximum  of  16.8.  The  final  reduced  order  model 
is  given  by  am  =  —  .7485,  bm  =  .0944,  and  cm  =  .9.  In  the  unstable  case,  the  intial  values 
were  am  =  1,  bm  =  1,  and  cm  =  1  and  the  final  values  are  given  by  am  =  .9505,  bm  =  .0812, 
and  cm  =  1.1119.  In  this  case,  the  value  of  A  increased  from  2.8  to  a  final  maximum  of 
17.4. 

The  time  responses  to  a  step  input  are  given  in  Figs.  4  and  5  (p.  37).  A  comparison  of 
the  time  responses  shows  excellent  correlation  between  the  plant  and  model  outputs.  The 
frequency  responses  are  also  well  matched,  at  least  up  to  10  rad /sec.  These  can  be  seen  in 
Figs.  6  and  7  (p.  38).  Divergence  in  the  frequency  responses  is  to  be  expected,  since  no 
low  order  model  can  match  the  frequency  response  of  the  high  order  plant  at  sufficiently 
high  frequencies. 

Example  2.  Aircraft  with  structural  modes 

In  this  example,  an  eighth  order  plant  will  be  reduced  into  a  fourth  order  system.  The 
plant  is  the  longitudinal  system  of  the  Advanced  Supersonic  Transport  (AST)  along  with 
the  two  lowest  frequency  structural  modes  [11].  The  structural  modes  are  the  first  and 
second  fuselage  bending  modes.  The  system  is  of  the  form 


xp  —  Apxp  "i"  xp(0)  0, 

(61) 

yp  —  Cpxp, 

(62) 

xp  =  (v  a  6  q  x\  ii  x-i  £2)* 

(63) 

u  =  (6t  6t  6C  6a)* 

(64) 

where  the  matrices  Ap,Bp,  and  Cp  are  given  in  Table  1.  In  (63),  the  variables  on  the 
right  side  are,  respectively,  perturbed  speed,  angle  of  attack,  pitch  angle,  pitch  rate,  first 
fuselage  bending  mode,  its  derivative,  second  fuselage  bending  mode,  and  its  derivative. 
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The  quantities  on  the  right  side  of  (64)  are  the  control  inputs  from  the  elevator,  throttle, 
canard,  and  elevon,  respectively.  The  flight  condition  is  supersonic  cruise  at  Mach  2.5. 
The  units  for  the  airspeed  are  ft /sec,  and  the  angles  and  the  control  surface  deflections  are 
in  degrees. 

We  attempted  to  reduce  this  to  a  fourth  order  system  given  by 

im  =  Amxm  +  Bmu,  xm(0)  =  0,  (65) 


J/m  —  CrnXm , 


xm  =  (v  a  9  q)* . 


Table  1  Plant  matrices  of  the  AST  longitudinal  system 


-0.0127  -0.0136 
-0.0969  -0.4010 
0.0000  0.0000 
-0.2290  1.7260 

0.0000  0.0000 
0.0000  0.1204 

0.0000  0.0000 
0.0000  0.1473 


-0.0360  0.0000 

0.0000  0.9610 

0.0000  1.0000 
0.0000  -0.7220 

0.0000  0.0000 
0.0000  0.0496 

0.0000  0.0000 
0.0000  0.3010 


0.0000  0.0000 
19.5900  -0.1185 

0.0000  0.0000 
-12.0200  -0.3420 
0.0000  1.0000 
-44.0000  -1.2740 
0.0000  0.0000 
-7.4900  -0.1257 


0.0000  0.0000 
-9.2000  -0.1326 

0.0000  0.0000 
1.8420  0.8810 

0.0000  0.0000 
-4.0300  -0.5080 

0.0000  1.0000 
-21.7000  -0.8030 


J5,= 


C„  = 


/  0.0000  0.0194 

-0.0215  0.0000 
0.0000  0.0000 
-1.0970  0.0000 
0.0000  0.0000 
-0.6400  0.0000 
0.0000  0.0000 
V -1.8820  0.0000 

(100000 
0  1  0  0  0  0 

0  0  1  0  0  0 

0  0  0  1  0  0 


0.0000 

0.0000 

-0.0040 

-1.7860 

0.0000 

0.0000 

0.3660 

-0.0569 

0.0000 

0.0000 

0.1625 

-0.0370 

0.0000 

0.0000 

0.4720 

-0.0145 

0  0\ 

0  0  | 

0  0  I 

0  0/ 
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Note  that  in  this  case,  A  is  a  function  of  48  independent  variables  and  the  matrix  F 
in  (54)  is  of  dimension  24  x  24.  The  exponentials  in  (55)  and  (56)  were  evaluated  using 
the  built-in  matrix  exponential  routine  of  PC-MATLAB.  We  observed  that  because  of  the 
large  numbers  involved,  it  is  best  not  to  invert  the  matrix  in  (55)  to  get  the  matrix  in  (56), 
but  to  compute  it  directly  using  the  built-in  routine.  We  ran  the  multi-start  hill  climbing 
algorithm  with  three  different  starting  points,  viz.,  with  the  truncated  system  matrices,  and 
with  the  reduced  order  matrices  from  [11],  which  were  obtained  by  balancing  and  spectral 
decomposition.  The  final  time  T  was  taken  to  be  5  sec,  the  weighting  matrix  Q  was  the 
identity  matrix,  and  the  weighting  matrix  R  was  selected  as  the  diagonal  matrix  with  all 
diagonal  entries  equal  to  0.001.  Although  all  the  three  runs  yielded  local  maxima  for  A,  we 
obtained  the  best  value  of  A  with  the  initial  matrices  obtained  by  spectral  decomposition. 
In  this  case,  the  value  of  A  increased  from  2.6  to  36.3.  The  resuting  reduced  order  model 
is  characterized  by 


A 


m 


B 


m 


c 


m 


/ -0.0112 

-0.0023 

-0.0258 

-0.0003 

-0.1248 

-0.4222 

0.0139 

0.8712 

-0.2817 

0.0142 

0.0152 

1.0198 

V  -0.2862 

1.7511 

-0.0032 

-0.6893 

/  0.0027 

0.0215 

0.0013 

-0.0005  \ 

0.5117 

0.0022 

-0.1354 

-1.8368  | 

-0.1422 

0.0107 

0.0449 

0.0035  1 

\ -1.0341 

0.0054 

0.3664 

-0.0521  / 

/  1.0427 

0.0099 

0.0035  0.0103  \ 

-0.1416 

1.0082 

0.0028  - 

0.0059 

-0.0047 

0.0068 

0.9932  - 

0.0006  ' 

\ -0.1852 

0.0197 

0.0128  1.0255  / 

(68) 


(69) 


(70) 


The  original  unaugmented  plant  has  the  short  period  eigenvalues  at  0.6687  (unstable), 
and  —1.7755  (stable),  the  phugoid  eigenvalues  at  —0.0151  ±*0.0886,  the  first  fuselage  bend¬ 
ing  mode  eigenvalues  at  —0.7257  ±  *6.7017,  and  the  second  fuselage  bending  mode  eigen¬ 
values  at  —0.3122  ±  *4.4484.  The  eigenvalues  of  Am  are  given  by  —0.0046,  —0.0232, 0.7113, 
and  —1.7910.  The  correlation  between  the  plant  and  model  outputs  is  excellent  in  the 
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chosen  interval,  and  this  can  be  seen  from  Figs.  8-11  (p.  39,40),  where  the  responses  to  an 
elevator  step  input  are  plotted. 

Even  though  this  example  has  weakly  coupled  flexible  and  rigid  body  modes,  it  re¬ 
vealed  another  interesting  feature.  A  comparison  of  the  frequency  responses  in  [11]  using 
spectral  decomposition  and  balancing  demonstrated  the  superiority  of  the  spectral  decom¬ 
position  method  in  this  case.  As  a  comparison,  in  Table  2,  we  list  the  initial  and  final 
values  of  A  with  starting  points  obtained  from  the  truncated  system  matrices,  balancing, 
and  spectral  decomposition.  It  can  be  observed  from  Table  2  that  spectral  decomposition 
gives  the  best  initial  and  final  values  for  A. 


Table  2  Values  of  A  with  different  starting  points 


Initial  A 

Maximum  A 

Truncation 

0.1 

2.0 

Balancing 

0.3 

3.3 

Spectral  decomposition 

2.6 

36.3 

While  the  starting  point  obtained  from  the  spectral  decomposition  method  gave  the 
best  value  for  the  correlation  measure,  an  improvement  in  the  value  of  A  can  be  seen  for  all 
the  starting  points.  This  suggests  smother  use  of  our  method.  Using  our  algorithm,  we  can 
fine-tune  the  reduced  order  models  obtained  by  some  other  model  reduction  methods.  We 
are  currently  developing  an  optimization  algorithm  utilizing  certain  characteristics  of  A  as  a 
function  of  the  reduced  order  model  parameters.  It  is  hoped  that  these  characteristics  will 
allow  us  to  overcome  the  problems  associated  with  the  large  number  of  local  maxima  and 
the  sharp  rise  in  the  value  of  A  near  the  global  maximum.  Global  optimization  algorithms 
currently  being  considered  include  stochastic  search  methods  [9],  methods  of  global  increase 
such  as  simulated  annealing  [9,12],  and  methods  of  improvement  of  local  maxima  such  as 
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tunneling  [9]. 

As  is  mentioned  at  the  beginning  of  this  section,  a  vast  amount  of  computation  time 
was  used  up  in  the  evaluation  of  A  for  given  reduced  order  model  matrices.  An  important 
topic  for  further  research  is  to  devise  an  alternate  method  for  the  efficient  evaluation  of  A. 
Also,  our  evaluation  of  A  is  only  accurate  to  within  the  specified  incremental  step  size  of 
A  and  this  may  lead  to  premature  termination  of  the  optimization  routine  that  seeks  to 
maximize  A. 


5.  Conclusions 

In  this  chapter,  we  presented  a  technique  for  reduced  order  modeling  with  a  modi¬ 
fied  Hoo  optimality  criterion.  A  characterization  for  the  determination  of  the  correlation 
between  plant  and  model  outputs  is  given.  Also  the  problem  of  robust  model  reduc¬ 
tion  is  addressed  and  an  expression  for  the  variation  of  correlation  with  plant  parameter 
variations  is  derived.  Nonlinear  programming  algorithms  were  utilized  to  reduce  the  lon¬ 
gitudinal  flexible-body  model  of  an  Advanced  Supersonic  TVansport.  Further  work  needs 
to  be  done  to  devise  a  suitable  global  optimization  algorithm. 
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Fig.  4  Step  Responses  in  the  Simple  Stable  Case 


Fig.  5  Step  Responses  in  the  Simple  Unstable  Case 
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Fig.  8  Angle  of  Attack  due  to  Application  of  an  Elevator  Step 


Fig.  9  Pitch  Attitude  due  to  Application  of  an  Elevator  Step 
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Fig.  10  Pitch  Rate  due  to  Application  of  an  Elevator  Step 


Fig.  1 1  Velocity  due  to  Application  of  an  Elevator  Step 
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